#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <vector>
#include "es-bvs-bl-common.h"

//#define DEBUG

//
//  benchmarking program
//
int main( int argc, char **argv )
{    
    if( find_option( argc, argv, "-h" ) >= 0 )
    {
        printf( "Options:\n" );
        printf( "-h to see this help\n" );
        printf( "-n <int> to set the number of particles\n" );
        printf( "-o <filename> to specify the output file name\n" );
        return 0;
    }
    
    int n = read_int( argc, argv, "-n", 1000 );

    char *savename = read_string( argc, argv, "-o", NULL );
    
    FILE *fsave = savename ? fopen( savename, "w" ) : NULL;
    particle_t *particles = (particle_t*) malloc( n * sizeof(particle_t) );
    double grid_size;
    set_size( n, grid_size );
    
    init_particles( n, grid_size, particles );
  
    const double cell_size = MIN_CELL_SIZE;
    const int bins_dim = (int)((grid_size+cell_size)/cell_size);

#ifdef DEBUG
    printf("grid_size = %lf\n",grid_size);
    printf("array of bins is = %d by %d\n",bins_dim, bins_dim);
#endif

    std::vector< particle_t* >* bins[bins_dim*bins_dim];
    for ( int i = 0; i < bins_dim*bins_dim; i++ )
    {
	bins[i] = new std::vector< particle_t* >();
    }
    
    int interact_neighbors[9] = {(-bins_dim-1), (-bins_dim), (-bins_dim+1),
				 -1, 0, 1,
				 (bins_dim-1), (bins_dim), (bins_dim+1)};
    
    //
    //  simulate a number of time steps
    //
    double simulation_time = read_timer( );
    int bin_x, bin_y, bin_num;
    std::vector< particle_t* > *inter_bin;
    std::vector< particle_t* >::iterator main_part, inter_part;
    for( int step = 0; step < NSTEPS; step++ )
    {
	//
	//  bin particles
	//
        for( int i = 0; i < n; i++ )
        {
	    bin_x = (int)(particles[i].x/cell_size);
	    bin_y = (int)(particles[i].y/cell_size);
	    bin_num = bin_y*bins_dim+bin_x;
	    bins[bin_num]->push_back(&particles[i]);
	}

        //
        //  compute forces
        //
        for( int i = 0; i < bins_dim*bins_dim; i++ )
        {
	    for ( main_part=bins[i]->begin(); 
		  main_part<bins[i]->end(); 
		  main_part++ )
	    {
		(*main_part)->ax = 0;
		(*main_part)->ay = 0;
		for ( int bin_to = 0; bin_to < 9; bin_to++ )
		{
		    // Edgar - if bin is on edge this will attempt redundant bin interactions
		    if (i + interact_neighbors[bin_to] < 0 ||
			i + interact_neighbors[bin_to] >= bins_dim*bins_dim )
		      continue;
		    inter_bin = bins[i + interact_neighbors[bin_to]];
		    for ( inter_part=inter_bin->begin(); 
			  inter_part<inter_bin->end(); 
			  inter_part++ )
		    {
			apply_force( **main_part, **inter_part );
		    }
		}
	    }
        }
        
        //
        //  move particles
        //
        for( int i = 0; i < n; i++ ) 
            move( particles[i], grid_size );
        
        //
        //  save if necessary
        //
        if( fsave && (step%SAVEFREQ) == 0 )
            save( fsave, n, particles, grid_size );
        
        //
        //  clean out the damn bins
        //
	for( int i = 0; i < bins_dim*bins_dim; i++ )
        {
	    bins[i]->clear();
	}
    }
    simulation_time = read_timer( ) - simulation_time;
    
    printf( "n = %d, simulation time = %g seconds\n", n, simulation_time );
    
    free( particles );
    if( fsave )
        fclose( fsave );
    
    return 0;
}
